NASA-CR-199035 


CTJ-CSSC-90-08 


April 1990 


V- 






CENTER FOR SPACE STRUCTURES AND CONTROLS 




A COMPUTATIONAL PROCEDURE 
FOR MULTIBODY SYSTEMS 
INCLUDING FLEXIBLE BEAM 
DYNAMICS 


(NASA-CR-199035) A COMPUTATIONAL N95-31730 

PROCEDURE FOR MULTI80DY SYSTEMS 
INCLUDING FLEXIBLE BEAM DYNAMICS 

(Colorado Univ.) 42 p Unclas 


G3/39 0060379 


by 

J. D. Downer, K. C. Park, and J. C. Chiou 


COLLEGE OF ENGINEERING 
UNIVERSITY OF COLORADO 
CAMPUS BOX 429 
BOULDER, COLORADO 80309 





A COMPUTATIONAL PROCEDURE FOR 
MULTIBODY SYSTEMS INCLUDING 
FLEXIBLE BEAM DYNAMICS 


J. D. Downer, K. C. Park, and J. C. Chiou 


Department of Aerospace Engineering Sciences and 
Center for Space Structures and Controls 
University of Colorado, Campus Box 429 
Boulder, Colorado 80309 


April 1990 


Report No. CU-CSSC-90-08 







A Computational Procedure for Multibody 
Systems Including Flexible Beam Dynamics 
J. D. Downer, K. C. Park, and J. C. Chiou 
Department of Aerospace Engineering Sciences 
and Center for Space Structures and Controls 
University of Colorado at Boulder 
Boulder, CO 80309-0429, USA 


Abstract 

A computational procedure suitable for the solution of equations of motions for flexible 
multibody systems has been developed. The flexible beams are modeled using a fully non- 
linear theory which accounts for both finite rotations and large deformations. The present 
formulation incorporates physical measures of conjugate Cauchy stress and covariant strain 
increments. As a consequence, the beam model can easily be interfaced with real-time strain 
measurements and feedback control systems. A distinct feature of the present work is the com- 
putational preservation of total energy for undamped systems; this is obtained via an objective 
strain increment/stress update procedure combined with an energy-conserving time integra- 
tion algorithm which contains an accurate update of angular orientations. The procedure is 
demonstrated via several example problems. 


1. Introduction 

The simulation of flexible multibody systems is becoming an increasingly important 
tool for the design and operation of many engineering applications. Typical examples of 
such systems include deployable space structures, high precision machine dynamics and 
robotics, and other problems containing controlled positioning of structural components. 
The components of these articulated structures typically undergo large relative displace- 
ments and rotations in order to carry out the intended operations. To perform the desired 
kinematic motions, various types of mechanical joints are introduced to constrain the rel- 
ative motion between the various components. New technology needs of both the space 
and robotics industries have increased the demand for accurate numerical simulations of 
the effect of component flexibility on the performance of multibody systems. A significant 
coupling between the gross structural motion and the elastic deformation can be expe- 
rienced by typical applications in which lightweight structures with higher flexibility are 
required to operate with greater positioning accuracy and at higher speeds. To capture 
this phenomenon, a realistic mathematical model of the structural component that can 
readily be incorporated into a general multibody dynamics methodology is necessary. 

Two basic approaches, the floating frame approach and the nonlinear continuum ap- 
proach, exist for the modeling of flexible components within a general multibody system. 


1 



.The floating frame approach introduces a moving reference frame to follow some overall 
mean rigid body motion of the beam; the elastic deformation of the beam is then de- 
scribed relative to this moving reference 1-6 . With this approach, the classical multi-rigid 
body analysis was extended to include structural flexibility by superposing existing linear 
deformation descriptions onto the rigid motions of the floating reference frame 7 ’ 8 . The 
definition of such a mean axis system and the corresponding deformation modes within 
the general context of the finite element method has been presented 9-11 . To minimize the 
number of elastic coordinates, coordinate transformations from the physical elastic coor- 
dinates to modal coordinates were performed within the multibody dynamics context 12 , 
and static correction modes were used in conjunction with the normal modes of vibration 
to account for reaction forces and torques transmitted to the components through joint 
connections 13,14 . An alternative choice of a floating reference frame for finite element appli- 
cations, termed the convected coordinate system, was introduced as a simple separation of 
the rigid body motion and the structural deformation for a given finite element 15-18 . All of 
these studies, however, are limited by the assumption of linear deformation theory which is 
inadequate to capture certain nonlinear phenomena. Nonlinear deformation theories must 
be taken into account for such instances as the geometric stiffening of a spinning beam 19,20 
in which the structural component experiences a centrifugal force as well as applications 
in which the components necessarily have low mass and very high flexibility. Extensions 
of the original approach to model the nonlinear effects include the substructuring tech- 
nique in which the component is further partitioned into substructures each with a local 
reference frame where normal vibration and static correction modes can then be used to 
model the deformation 21 , and the finite element incorporation of a nonlinear Green strain 
measure 22,23 . The resulting equations of motion of the floating frame approach, written in 
terms of a set of reference coordinates and a set of relative elastic coordinates, inherently 
contain a complex coupling of the gross motion and the elastic deformation modes. 

Recently, a fully nonlinear continuum approach to describe the dynamics of the flexible 
beam has been pursued 24-28 . Through the use of finite-deformation rod theories 29-32 , the 
approach is capable of directly accounting for both finite rotation kinematics and large 
deformations of the beam component. Since the motion due to rigid rotations of the beam 
is not distinguished from that due to deformations, the need for a floating reference frame 
is completely obviated and the component inertia is identical in form to that of a rigid 
body. The inherent nonlinear coupling between the gross body motion and the elastic 
deformation is transferred to the stiffness part of the equations of motion. The key to 
the successful adoption of this approach is to develop a computational procedure for the 
nonlinear internal force term that preserves rigid body motions. 

The aim of this paper is to incorporate the nonlinear continuum formulation of the 
spatial beam motion into a general multibody dynamics software methodology. The present 
formulation employs a convected coordinate representation of physical Cauchy stresses 
and corresponding set of physical strains. This representation naturally lends itself to the 
“software in the real-time experiment” loop as sensors measure only physical quantities. 
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Another advantage of the formulation is that the degrees of freedom of the beam component 
embody both the rigid and flexible deformation motions. The task for incorporating the 
multibody system constraints becomes straightforward, and the equations of motion for an 
arbitrary configuration of flexible beams and rigid bodies can automatically be generated 
in terms of an identical set of physical coordinates. Numerical solution procedures for 
the integration of spatial kinematic systems can then be directly applied to these physical 
coordinates. Such a universal treatment is not applicable within the context of the floating 
frame approach as the reference and elastic coordinate definitions are of highly different 
character. 

The rest of the paper will be organized as follows. Section 2 will detail the kinematic 
description of the continuum beam in which the total motion is referred directly to the 
inertial reference frame. The principle of virtual work of a continuum as specialized to 
the spatial motion of the beam component is detailed in Section 3. The subsequent finite 
element discretization of the beam component and overall multibody system equations are 
then presented. Section 4 will summarize the staggered procedure for the integration of 
multibody dynamic systems. The virtual work expression is used to derive the method 
of computation of the internal force, and Section 5 will address this algorithmic treat- 
ment of the nonlinear stiffness operator. Section 6 will present some example problems 
demonstrating the software capabilities. 


2 . Beam Kinematics 

The present formulation adopts an inertial reference frame for describing the trans- 
lational motions and a body-fixed frame for the rotational motions. The consequence of 
this description is that the translational and rotational variables embody information due 
to both rigid rotations and deformations of the beam. The configuration of the beam, as 
shown in Figure 1, is completely characterized using a position vector locating the neutral 
axis of the beam from the inertial origin and a body-fixed frame representing the orienta- 
tion of the cross-section with respect to the inertial reference frame. The position vector 
r locating an arbitrary particle point on the beam is thus described as 

r = ( X + u ) T e -I- b ( 2 . 1 ) 

where “boldface” symbols represent three subscripted vectors and the normal type symbols 
represent three components of a given vector; e = { ei,e2, e3 } T represents the three 
orthogonal vectors defining the inertial reference frame; b = { bi,b2,b3 } r represents 
the body-fixed reference frame which is attached to and rotates with the beam cross section; 
X = { Xi, .Y 2 , X 3 } T represents the inertial components of the original neutral axis 

position; u = { Ui,u 2 ,«3 } T represents the inertial components of the subsequent total 
translational displacement of the neutral axis, and l T = { 0, £ 2 >^3 } are the body-fixed 
components of the distance from the beam neutral-axis to the material point located on 
the deformed beam cross-section. It is noted that the beam cross-section is allowed to 


3 



rotate such that it is not necessarily perpendicular to the neutral axis in order to model 
transverse shear deformations. Warping deformation of the cross-section is not taken into 
consideration. 


In order to derive the necessary time derivatives for the description of the large rotation 
dynamics, we employ the well known formula 33 : 

d d e d b 

di = dt = dt +UJX 

where oj is the angular velocity vector and the superscripts e and b indicate that the 
derivatives are to be those observed in the inertial (space) and body (rotating) system of 
axes respectively. The above is expressed in the matrix form to act on the body frame 
components of a given vector 

d d b 

Tt ~ Tt + " ’ (2 3) 

and the velocity and acceleration of the position vector (2.1) are 


dt 

d 2 u T 

~dt r 


e + i T u T b 




+ Q T uf T ) b 


Given the following relation between the b-basis and the e-basis 


b = R e 


where Ris a (3 x 3) orthogonal transformation matrix, the body frame components of the 
skew-symmetric angular velocity tensor ( u) T ) are 


~ T 

U 7 = 


A conjugate virtual rotation tensor is defined analogous to the above as 

Sa T = <5RR r , 

and the variation of the position vector (2.1) is given as 


<$r = Su T e + £ I 6a I h . 


T c~T\ 


The equations of motion as derived from the stated beam kinematic description will be 
discussed next. 
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3. Spatial Beam Equations of Motion 

The principle of virtual work, which is simply a ‘weak’ or variational form of Cauchy’s 
differential equations of motion for the equilibrium of a given set of particles of a continuum, 
is stated as 34 

J v Sn pfi dV + atj d £dV = J y Snft dV + J s SnU dS . (3.1) 

The cartesian coordinates x, represent the particle position after some deformation has 
taken place, <5r, a kinematically admissible virtual displacement, f; the acceleration, /, 
the external force per unit mass, and t{ the stress vector acting on a surface with outward 
normal components rij. Likewise, represents the cartesian components of the Cauchy 
stress tensor, and p is the mass density. The expression is tailored for the continuum 
beam by using the kinematic relations (2.1), (2.4), and (2.8) for the components x,, Sri, 
and ?i respectively. As well as providing the basis for a finite element approximation 
techniques, the variational formulation readily lends itself to the derivation of incremental 
strain-displacement relations as deduced from the derivatives of the virtual displacement 
components. The present formulation employs a physical stress measure defined as a force 
per unit deformed area and the conjugate physical strain increments based on the de- 
formed coordinates. As such, the formulation can be recast into a convected coordinate 
system moving with the beam, thus simplifying the stress and strain computational proce- 
dures. The practical advantages of such a formalism axe in real-time software simulation 
experiments as the computed physical quantities correspond to the actual stress/strain 
measurements of the sensors located and operating on the deformed structure. 

For notational convenience and subsequent finite element discretization, the principle 
of virtual work is expressed in the following operator form: 

8F 1 + 8F S = 8F e + 8F t (3.2) 

where the inertia operator 8F 1 , internal force operator 8F S , external force operator 8F e . 
and traction operator 8F T are identified from (3.1). Explicit expressions for the various 
operators incorporating the large rotation beam kinematics are derived in Sections 3.1 to 
3.3. The finite element discretizations are given in Section 3.4, and the incorporation of 
the beam formulation into the multibody dynamics framework is discussed in Section 3.5. 


3.1 Spatial Beam Inertia Operator 


The inertia operator was defined from (3.1) as 

8F 


1 = IP 8riri dV = / p 8r • r dV 

Jv Jv 


(3.3) 
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from which an expression can be derived directly from the kinematic equations (2.4) and 
(2.8). If the origin of the body-fixed basis is located at the centroid of the cross-section, 
the following simple expression results for SF 1 : 


6F 1 


-I. 


{ Su T Sa T } < 


pA 


di * 


ds 


T 

J dt 


(3.4) 


-f- Co J u> 


where 



p IF dA 


J 


represents the inertia tensor of the beam cross-section and ds represents the remaining 
integration to be performed over a beam length parameter. The translational inertia is 
completely decoupled from the rotary inertia and is of the same form as that seen in rigid 
body dynamics. This is due to the dual choice of the translational displacements measured 
in the inertial basis and the angular velocity measured in the body-fixed basis located at 
the center of mass of the cross-section. 


3.2 Spatial Beam Internal Force Operator 

The internal force operator was defined in (3.1) as 

4F$ = jy ST ^ ^ (3 5) 

identifying as conjugate quantities the virtual displacement gradient and the Cauchy stress 
tensor. This form of the internal force along with the beam kinematic description will be 
used to deduce a set of virtual strain-displacement relations that are invariant to rigid body 
motions. The corresponding conjugate stress tensor will be obtained from an objective 
incremental procedure that relates incremental strains obtained from the virtual strain 
tensor to Cauchy stress increments. Thus the internal force term will be derived completely 
from the original definition of the beam kinematics without making an a priori definition 
of the existing strains or stresses. 

A physically appealing decomposition of the stress and virtual strain tensors into an 
alternative beam reference frame which lies tangent to the deformed neutral axis is intro- 
duced to provide conceptual simplifications in the derivation and subsequent computations. 
Certain stress states referenced to this convected frame are kinematically required to van- 
ish in a beam formulation. When applied to the convected frame stress components, this 
choice also leads the task of stress update computations to a simple additive procedure. 
To this end, we introduce a convected reference frame, denoted by a, which is related to 
the inertial reference frame e by 

a = T e , a = { a!,a 2 ,a 3 } T . (3.6) 
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For implementation purposes within the context of the finite element method, the con- 
vected frame will be constant on the element level and thus is similar in concept to that 
introduced by Belytschko et al. 15,16 . It is noted that this reference frame does not coincide 
with the body frame b attached to the cross-section. The relative difference between these 
two reference frames is represented by the rotation matrix S which models the effects of 
transverse shear and torsion deformation as 


b = Sa, R = ST , 

and the latter interdependence between the rotation matrices is established. 


(3.7) 


The internal force operator, originally characterized by the inertial frame components 
of the Cauchy stress tensor ( ) and conjugate virtual displacement gradient, will equiv- 

alently be expressed in terms of the convected frame components of the stress tensor ( a*- ) 
and a corresponding convected virtual displacement gradient as 


8F i 




98 r, 

dV 


-L 


T „« AV 

*mt a mk dv 


(3.8) 


The symmetric portion of the transformed deformation gradient is used to define the virtual 

d8ri 


strain tensor 8e^ nk as 


6e a = — ( T 

mk ~ 2 ' m ' d£ k 


+ Tfci 


din 


) 


(3.9) 


which is an objective tensor invariant to arbitrary rigid body motions. The internal force, 
written in terms of the convected frame tensors, will be expressed in vector format as 


6F i 


-L 


8e a mk * a mk dV 


f ( 1 

= J { % } S tejr, > 

V ( Se K ) 


dV 


(3.10) 


where the notation 


£* — { £2* £3 } — { £?* 7 >C } 1 ( )ij — ( )ij + ( )ji 

denotes the coordinates of the convected reference frame and the engineering shear strain 
definitions respectively. The rest of the convected frame strain components 

are identically equal to zero due to the original assumptions of the beam kinematics. 

A set of virtual strain-displacement relations can be derived from the expressions (2.S) 
and (3.9). The final result is expressed as 



= £7 + r 6 k 


(3.11) 
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where 


d 68 rp 

6k = , 6/3 = S T 6a (3.12) 

and is comparable to that of Reissner 31 . In the above expressions, 67 represents the 
membrane and two transverse shear strains, 6k the torsion and two bending strains, and 
6/3 the virtual rotations of the cross-section referred to the convected frame. 


. _ d6u . 

= T ir + < 


0 

-60 3 
6(3 2 


In an analogous manner the total stress state is expressed as 

' 

> 

to be obtained from a separate stress update procedure. A substitution of (3.11) and (3.13) 
into (3.10), and a spatial integration over a symmetric cross-sectional area results in the 
following expression for the internal force 


= a. 


+ P 


(3.13) 


6F S = J^{ 6 7 t N y + 6 k t M k } d£ (3.14) 

where Ny represent the axial and transverse shear forces per unit length, and M K represent 
the torsional and bending moments per unit length as given by 


Ny = f <7 dA , M k = f i 1 " t 7 dA . (3.15) 

J A dA 

To be consistent with the inertia operator derived in (3.4), the above is written as 

6 F S = J { Su * Sa T } [ B } T | ) d( (3.16) 

which involves a transformation back to the body frame components of the virtual rotations 
and also an identification of the desired incremental strain-displacement matrix B. To 
effect the change of the body reference frame of the cross-section orientation in space with 
respect to the constant convected reference frame, we invoke the following relations: 


d a 6/3 

dt 


gT d a 6& _ gT ^ d b 6ot 


d( 


dt 


+ ks 6a ) 


~T d*S T 

= w s 


(3.17) 
(3. IS) 
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which axe completely analogous to those relating changes in the time derivative given in 
(2.3) and (2.6). The strain operator [ B ] of (3.16) is then recognized as 


r 

J- at 

i i 

S T 

*1 


' 0 

0 

0 ' 

a? 



. 

9 H = 

0 

0 

-1 

0 

S T ( K S 

+ I 

*)J 


0 

1 

0 


(3.19) 


It remains to provide a procedure for updating <r 7 and a K in order to compute iV 7 and 
M k . For this purpose, we employ the following rate-type law that relates the instantaneous 
rate of stress to the instantaneous rate of deformation: 

°kl ^ Cklmp £ a mp (3.20) 

where Ckimp represents the material response tensor, and and imp represent the 

convected frame stress and strain rates, respectively. This approximate constitutive law can 
be derived by transforming the Truesdell rate equation 35 , which is an objective equation 
based on inertial components of Cauchy stresses and the velocity gradient tensor, to the 
convected basis. This equation is then integrated in time as 


<7 


a 

kl 


n + l 




/*t n + l 

a " 
kl 

+ 

/ c klmp £ mp dt 



Jt n 

a n 
kl 

+ 

Ckimp Ae mp 


(3.21) 


to define the stress update procedure. The evaluation of the strain increments to 

be defined from the virtual strains (3.12), will be detailed in Section 5. 


3.3 Spatial Beam External Force and Traction Operator 

The external force operator defined in (3.1) as 

6F e = f Sri fi dV 
Jv 

has the final resultant form 

6F E = j { 6u T 6a T ) { f '„ } d( (3.22) 

where f e represents the inertial components of a force per unit length acting on the beam 
neutral axis and f b represents the body-fixed components of a moment per unit length 
acting on the beam cross-section. The traction operator defined as 

6F t = j Sri t { dS (3.23) 

Js 
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acts on the exterior surfaces of the beam as natural boundary conditions. 


3.4 Finite Element Discretization 

The variational form of the partial differential equations representing the spatial dy- 
namics of a continuous beam presented in the preceding sections provide a basis for the 
finite element method to be used as a spatial discretization procedure 36 . In the present 
study, we restrict ourselves to the use of linear shape functions to approximate the dis- 
placement field along the beam, viz., 


npe 

u — N i uj (3.24) 

/=i 

where Nj denotes the spatial linear shape functions, u/ represents the degrees of freedom 
at the element nodes, and npe denotes the number of nodes per element. The element 
inertia operator, from (3.4), is written as 


npe npe 


d 2 u i 


d b uiu' 


6F 1 — { Suf pA M E IK ^ t 2 + Saf pJ T M E 1K —— } 


/=1 R'=i 


dt 


npe 

+ dE i 

1=1 


(3.25) 


where 


m e 


IK 


= / Nj Nk d£ , 


D E (u)t 


= J ( uJu ) 


I d£ 


represent the element mass matrix and nonlinear angular acceleration vector. The former 
will be evaluated as a standard lumped mass matrix for the computational efficiency of 
explicit integration techniques to be described in Section 4, and the latter will be evalu- 
ated from an average of the element nodal angular velocities. The element internal force 
operator, from (3.16), is written as 


SF s = 


" Pe r ( M } E npe ( Ce 'i E 

E {"£} = E {*<•; *“'>{ 5 ;} ( 3 - 2 « 


where the evaluation of the element strain operator 


[*'] = l 


*1 Nj S J 

0 SjiksNr + 2£)\ 


dt 


(3.27) 
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and the resultant element stresses N y and M K , as defined in (3.19) and (3.15) respectively, 
will be presented in detail in Section 5. The element external force operator, from (3.22), 
is written as 


8F 


E 


= £ {**/ *>/) {/j} . fi* = l N ‘f' ,b 


di 


(3.28) 


and the traction operator is implemented as boundary conditions on the nodes. The 
equations of motion in terms of nodal degrees of freedom ( 8u d , 8ot d ) for the entire beam are 
obtained from an assembly of the above element operators. For the unconstrained beam, 
these nodal virtual displacements and rotations are arbitrary independent variations, and 
the discrete equations of motion are written as 


M d 0 
0 J d 


u d 


“>d 


+ 


i D d (u) J 


+ 



U (3 - 29) 


where M d , J d represent the assembled mass and inertia matrices, and D d (u>), S d , f d 
represent the assembled nonlinear acceleration, internal force, and external force vectors 
respectively. 


3.5 Extension to Multibody Dynamics 

The present formulation of spatial beam dynamics as given by (3.29) can readily be 
incorporated into a general multibody dynamics methodology. The degrees of freedom of a 
rigid body, namely the inertially-based translational position of the center of mass and the 
rotational orientation of the body reference frame, coincide with the degrees of freedom 
of the nodal coordinates of the present beam components. Thus the equations of motion 
(3.29) can be specialized to represent a rigid body system by setting the internal force S d 
equal to zero. 

It remains to augment both the holonomic and nonholonomic constraint conditions 
modeling the contacts among the various bodies to the equations of motion. For this pur- 
pose, the Lagrange multiplier technique is used to couple the algebraic constraint equations 
with the differential equations of motion of the generalized coordinates by augmenting the 
virtual work of the unconstrained system (3.2) with the virtual work requii’ed to enforce 
the constraints. Given a set of equations representing holonomic constraint conditions 
between the displacement coordinates as 

$ h ( u,t ) — 0 , 8$ h ~ — r— — ■ 8u = Bh 8u = 0 (3.30) 

au 

and a set representing nonholonomic constraint conditions between the virtual displace- 
ments and rotations as 

8$n ( u, 8u, R, 8a ) = = 0 , (3.31) 
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the virtual work expression (3.2) of the unconstrained system is modified to account for 
the constraint via Lagrange’s multipliers A as 37 

SF i + SF ^ -f- A h'S*&h + A = SF E + SF^ 


The virtual displacements and rotations of the generalized coordinates can now be treated 
as arbitrary independent variations in the modified virtual work expression. The equations 
of motion for constrained flexible multibody systems with respect to a set of generalized 
coordinates ( ii,u> ) denoting both the nodal coordinates of the flexible members and the 
physical coordinates of the rigid bodies can be expressed as 


where 


M 0 f u 1 

0 J \ u) J 


+ B r A 



f Qu 1 _ / f e -S e (u,q) \ 

\ QwJ \/ 6 -Z>( W )-S 6 («,?)J 


B 


Bh 

Bn ’ 


A 


(3.32) 

{£} 


in which D{ui) represents the nonlinear acceleration, S the internal force vector, / the 
external force vector, and B T A the constraint force vector. As an additional number of 
unknown Lagrange multipliers A for each constraint condition have been introduced along 
with the generalized coordinates for each degree of freedom, the above system of equations 
must be augmented with the constraint equations themselves to achieve a determined 
system of equations. 


4. Time Integration Techniques for Constrained Systems 

The present methodology to formulate the equations of motion of an arbitrary assem- 
blage of interconnected flexible beams and rigid bodies is readily adaptable for use with 
existing multibody dynamics solution techniques. The equations (3.32) model the beam 
components with degrees of freedom u and u j that embody both the rigid and flexible 
deformation motions. As such there is no need to solve separately generalized coordi- 
nates denoting the flexible motion from a reference set of coordinates denoting the rigid 
motion. In addition, as the nodal coordinates of the beam components are defined in 
the same kinematic manner as the physical coordinates of the rigid body components, no 
distinction need be made between the treatment of the flexible and rigid components of 
the multibody system other than the calculation of the internal force of the flexible mem- 
ber. Therefore, the salient feature of this type of formulation is that numerical solution 
procedures for the integration of spatial kinematic systems can be directly applied to the 
generalized coordinates of both the rigid and flexible components. 

A multibody dynamics solution procedure, originally demonstrated on rigid body sys- 
tems in previous studies 38-41 , is adopted for the above flexible multibody system equations 
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of motion. The key to the procedure is a staggered implementation of the separate gener- 
alized coordinate integrator and constraint force solver modules. An improved variation of 
the explicit central difference algorithm, described in Section 4.1, is used to integrate the 
translational displacements and the angular velocity of the system. An algorithm based on 
the Euler parameter representation of finite rotations, described in Section 4.2, is used to 
update the configuration orientation from the angular velocity. The computations of the 
Lagrange multipliers are then carried out in a separate routine, described in Section 4.3, 
which implicitly integrates a stabilized companion differential equation for the constraint 
forces in time. 


4.1 Explicit Generalized Coordinate Integrator 


The central difference explicit integration algorithm is written as 


d"+* = d n ~2 + hd n 

d n+1 = d n + hd n+ > (4.1) 

d n+1 = M~ l Q ( d n+l , d n+1 ) 

where the superscript n = 1,2,3, ••• designates the discrete time station t n = n h and 
h is the stepsize. Unlike in conventional structural dynamics, a straightforward application 
of (4.1) on the rotational equations 


J UJ + U) J U) — f u 


inherent in the multibody system equations of motion (3.32) leads to computational dif- 
ficulties. In order to compute a;”" 1 " 1 , it is necessary to have u; n+1 . However, due to the 
inherent nature of the algorithm, only u > n+ 2 is available. It was shown 41 that the naive 
approximation 


u 


n+l 


~ W n+ 2 


(4.2) 


results in a computationally unstable integration of the singular velocity u >. To correct 
this within the context of explicit computational sequences, an interlaced application of 
the central difference algorithm such that the displacements and velocities are advanced 
one- half time step at a time was proposed 40,41 . The algorithm advances the solution to 
the time station given the solutions of the two preceding time stations t n ~ i and t n 

as follows: 
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(«) 

lt n+ J 

= 

+ hu n 

(*) 

ti n+ 2 

= 

u n_ 2 + hu n 

(c) 


— 

a;" - * + hu n 

(d) 

q n+ i 

= 

q ( w n+ * ) 

(e) 

s n+ * 


S( u n+ s 9 n+ * ) 

(/) 


— 

D ( u n+ > ) , / n+ * = / ( i n+ * 

(5) 

Q n+ * 

= 

Q ( / n+ *,S n+ *,D n+ * ) 

(h) 

A n+ * 

= 

A ( A n ,Q n+ * ) 

(0 

u" + * 

= 

M" 1 ( - Bl A" + * ) 

U) 


— 

J” 1 ( - Bl A" + ± ) 


The evaluation of the generalized rotational parameters q to be obtained from the angular 
velocity, as represented by step ( d ), will be detailed in Section 4.2. The evaluation of the 
interned force 5 from the current configuration coordinates u and q, as represented by 
step (e), will be detailed in Section 5. The evaluation of the Lagrange multipliers A, as 
represented by step (h), will be detailed in Section 4.3. The algorithm proceeds to the next 
half time station < n+1 , now given the solutions at time stations t n and t n+ 2 , and thus the 
force and acceleration terms are evaluated twice each time step. The algorithm is initiated 
for time t * given initial conditions for time t° in the following manner: 


(k) 

u* 

= u° 

+ 

h ,. 0 

2 U 

(0 

u,* 

= u ; 0 

+ 

h . o 
2" 

(m) 

H 1 

= u° 

+ 

h 

(n) 

1 

U2 

II 

O | !-• 

u° 

+ u 1 ) 


from which steps ( d ) through (j ) can be performed. 

One last remark will be made on the angular velocity integration. The equations of 
motion were derived using body frame angular velocity components. The integration of 
these quantities shown in step (c) is not formally correct as the components at different 
time steps are defined with respect to different body-fixed frames. This concern can be 
eliminated by applying the central difference update to the inertial components of the 
angular velocity. Step (d) will then consist of an appropriate function of inertial angular 
velocity components. The integrated inertial angular velocities must be transformed to the 
moving reference frame before evaluating steps (/) and ( j ) since the equations of motion 
are written with respect to the body frame angular velocity description. The angular 
acceleration evaluated in step (j) must then be transformed back to inertial reference 
frame before being integrated again in step (c). 
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4.2 Rotational Parameter Integration 


The two-stage explicit integrator was applied to the translational displacement and 
velocity coordinates and the angular velocity coordinates. As the rotational orientation 
parameters axe not directly integrable from the angular velocity vector, a procedure must 
be developed to update the configuration orientation given the angular velocity. Any finite 
rotation can be uniquely expressed by a rotation angle 6 and an appropriate rotation axis 
n 42 . Two rotational parameterizations based on this description are the rotational vector 
( 0 ) and the Euler parameters ( g 0 , q ) defined respectively as 


0 = 9 n 


f 9o \ _ J cos | 1 

l q J l » sin | / 


(4.3) 


The three parameters of the rotational vector are independent, while the four Euler pa- 
rameters axe subject to the constraints 


9o + q = 1 

The rotation matrix is represented as a function of the Euler parameters as 


R = 


2 ( 9 o + 9 i )~ 1 2(9192 + 9093) 2(9193 — 9092) 

2 ( 9 i 92 ~ 9093) 2(90 + 92) — 1 2(9293 + 9091) 

2(9i 93 + 9092) 2(9293 — 9o9i) . 2(9 q + 93) ~ 1 


The body frame components of the angular velocity tensor defined in (2.6) as 


- T 
W 6 


= rr t = 


0 

U>3 



-CU3 

0 


, Ub = < u>2 j 

U>2 

-U>1 

0 

[ UJ 3 J 


has the Euler parameter representation 42 

- ■[ -• 


9o 

q 


q 

9oI-q 


]{« 


A similar expression for the inertial components of the angular velocity tensor 

uj = R T d>fR = R T R 

can be derived as 


0 } 


= 2 

1 We J 

. 


9o 

-q 


q 

9oI + q 




The above definitions can be inverted to yield the expressions 


It} ■ H 


0 

Ub 


-ul 


cot 


{:} - 


Ab { u>b ) 


{:} 


( 4 . 4 ) 


( 4 . 5 ) 


( 4 . 6 ) 


( 4 . 7 ) 


( 4 . 8 ) 


(4.9) 
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for the body frame components and 



1 

' 0 -uj ' 

(M 

= A* ( w e ) ■ 

f go 

l q J 

~ 2 

_ W e Cj t 

1 q 1 


l q 


for the inertial frame components. A general representation 


q = A(u>) q , * = { q } 


(4.10) 


(4.11) 


will be used to denote (4.9) or (4.10) given the angular velocity description. These in- 
verse expressions are derived from (4.6) and (4.8) by incorporating the derivative of the 
constraint equation (4.4) 

go go + q T q = 0 . (4.12) 


The configuration orientation is obtained from a numerical time discretization of the 
above Euler parameter - angular velocity representations. Among several possibilities, the 
approximation that satisfies the constraint condition (4.12) in the discrete sense is the 
following trapezoidal formula 

J ( 1 n+1 - «“ ) = A(«" + i) i ( ,“ +1 + «* ) (4.13) 

Due to the structure of A, the solution matrix can be analytically inverted such that the 
discrete orientation update 

,"+■ = I[/+ jAKli)] [/+ f A(u>-+*)] 4 - (4.14) 

where 

D = 1 + ~ ( w\ + + wf ) 

results. The final result is normalized to satisfy the constraint (4.4). The above equation 
is valid for either the body or inertial frame decomposition of the angular velocity as long 
as the corresponding form of A from (4.9) or (4.10) is used. The resulting update (4.14) 
involves only explicit computations and is readily incorporated into the two-stage explicit 
integration procedure. 


4.3 Constraint Force Solution Procedure 

A partitioned solution procedure has been employed to solve the generalized coor- 
dinates separately from the Lagrange multipliers. To effect a partitioned solution of the 
constraints, a stabilized companion differential equation for the constraint forces is formed 
by adopting the penalty procedure 38,39 . The penalty procedure uses the equations 

\ H = - $ H \ N = - $ N e — * 0 (4.15) 

e e v 
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as the basic constraint equations instead of (3.30) and (3.31) for the holonomic and the 
nonholonomic constraint conditions respectively. The penalty equations can be written in 
the general form, from (3.30) and (3.31), as 

* = \ Bi - i = {!} < 4 - 16 > 

The numerical solution to the above companion differential equation is obtained as follows. 
The constrained equations of motion (3.32) are integrated once from (3.20) using the 
implicit integration rule 

i"+i = i n + 8 z n+ ± , 8 = ^ 

£ 

as 

i n+ * = 8 M~ l ( Q n +* - B t A n+ * ) + z n (4.17) 

This expression is substituted into (4.16) to obtain the stabilized differential equation for 
the Lagrange multipliers 

c A n+ s + SB M~ 1 B t A n+ * = SB M~ l Q n+ * + B z n . (4.18) 

The same integration rule is applied to this equation to result in the discrete update 

(el + 8 2 B M~ l B T ) A n+ * = eA" + r" + * (4.19) 

r” + * = 6 2 B M~ 1 Q n+ ^ + 8 B z n . (4.20) 

The same procedure can also be derived with different integration rules. The update of 
the Lagrange multipliers, performed for each half time step, is easily adapted into the 
two-stage explicit integration procedure. 


5. Internal Force Computations 


The algorithmic treatment of the nonlinear stiffness operator is addressed in this 
section. The explicit generalized coordinate integrator of the previous section requires an 
evaluation of the internal force at a current time step t n + l from the coordinates of the 
beam configuration at that time. The internal force is first evaluated on the element level 
for all the finite elements comprising the flexible component from (3.26) as 



(5.1) 


after which these individual element computations are assembled to form the internal force 
of the discrete beam. The necessary computations to be described are the evaluations of 
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the discrete strain operator [ B E ] defined in ( 3 . 23 ) and the resultant stresses N y and 
M* respectively. 


The Timoshenko beam formulation in which the translational degrees of freedom are 
independent from the rotational degrees of freedom requires an approximation within the 
element such that these variables will be continuous across the element boundaries. Thus 
a two node finite element representing a linear interpolation of the translational and rota- 
tional variables is a sufficient discretization of the beam. To avoid the locking phenomenon, 
the interpolation of the rotational degrees of freedom associated with the transverse shear 
strain is underintegrated. After incorporating these concepts into ( 3 . 27 ), the resulting 
expression for the discrete strain operator is given by 



*uSf 

0 0 Sf ( - {I ) 


which acts on the virtual displacements and rotations 


2*1 Sj 1 

S 2 ( + 7 I ) 


( 5 . 2 ) 


{ 8 u\ 8 u 2 8 a.\ 8 a 2 } T 

where the subscripts refer to the element node number. The convected frame T matrix, 
body frame curvature tensor ks, and element neutral-axis length t axe constant quantities 
over the element domain, while the relative cross-section deformation S matrices are nodal 
quantities. The computation of these terms from the nodal displacement and rotation 
coordinates of the current configuration axe detailed in Section 5.1. 


A stress update procedure of the form 


N, 

M k 




( 5 . 3 ) 


is used to derive the resultant stresses of the current configuration at time t n+1 from the 
resultant stresses of the past configuration at time t n . The simple additive form of the 
procedure, which was derived from the numerical integration of a rate-type constitutive 
law, is due to the use of a convected frame stress and conjugate strain decomposition. The 
resultant stress increments AN-, and A M* axe obtained via 



' EA 

0 

0 ' 



'GJ 

0 

0 ' 

A iV 7 = 

0 

GA 

0 

A7 , 

AM* = 

0 

ei 2 

0 


0 

0 

GA 



0 

0 

1 


A set of strain increments A7 and Ak, which denote the change from time t n to t n+1 , are 
defined as a finite analogy to the infinitesimal virtual strains £7 and 8k derived in Section 
2 . A specific computational procedure designed for use with this incremental interpretation 
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of the continuum-based formulation such that the computed finite strain increments are 
invariant to arbitrary rigid body motions is discussed in Section 5.2. 


5.1 Computation of the Strain Operator 

The reference frames introduced in the formulation, namely the body frame b attached 
to the cross-section and the convected frame a tangent to the deformed neutral axis, are 
computed as follows. The Euler parameters representing the orientation of the beam cross- 
section at the finite element nodes are output from the generalized coordinate integrator 
at each time step. The rotation matrices R, , representing the b reference frame at each 
element node, are thus computed directly from the Euler parameter representation of a 
rotation matrix (4.5). This matrix contains rotational information of both that due to 
the rigid motion of the convected reference frame and the transverse shear and torsional 
deformations of the cross-section relative to the convected frame. 

The neutral axis of the finite element is defined as the straight line connecting the 
two element nodes, the tangent of which is trivial, and is directly calculated from the 
translational displacements output from the generalized coordinate integrator. Given this 
tangent ai, the a2 vector is defined as the cross product of aj with the b3 axis of Rj, 
and the remaining axis a3 defined to complete the right-hand coordinate system. The 
computed axis { ai , a2 , a3 }, as shown in Figure 2, define the rows of the T matrix. 
The rotation matrices S, , defined at each element node as the relative difference between 
the element convected frame and the nodal body frames, are thus 

Si = R, T t , i = 1,2 . (5.5) 

The procedure is an approximation applicable for moderate strains such that the S, matri- 
ces contain information solely due to transverse shear and torsional deformations 43 . The 
rotation matrices of the discrete strain operator (5.2) have thus been defined. 

The body frame components of the curvature tensor £5 defined in (3.18) as 

( Kx ' 

K = < K'2 > 

( *3 , 

are equivalent to 

rr d e R 'T 

*s = -gf K (5.6) 

as the convected frame T matrix is defined to be constant along the element domain where 
the differentiation is performed. This definition is completely analogous to the angular 
velocity tensor defined in (2.6) and motivates the use of an Euler parameter representation 
of the curvature completely analogous to the Euler parameter representation of angular 


r-T _ 

K S — 


- 

d( 


0 AC 3 — AC 2 

-“AC 3 0 K\ 

AC 2 — ACi 0 
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velocity (4.6) as a basis for the computation of the element curvature from the nodal 
rotational variables. The Euler parameter - curvature representation is 


?v O 

II 

to 

1 1 

9 o 

-q 

T 1 

q 

9oI - q . 

{1} = E(«.) 

dq 

ae 

(5.7) 

subject to the constraints 






% + q q 

= l 

dqo 

’ l( qo 

+ 

CdI 

£ 

II 

o 


(5.8) 

An approximation to be used in 

(5.7) 

that satisfies the constraint conditions in 

the discrete 

sense 

dq 1 

= 1 { 

- 9i 

) 5 Qa 

_ 2 ( 9l + 92 ) 


(5.9) 

If ( 9i + 92 ) 

II 


is evaluated using the Euler parameters of the element nodes output from the generalized 
coordinate integrator. It will be shown that this discrete computation is invariant to rigid 
rotations contained in the total nodal Euler parameters. 


The simple normalized average of the nodal Euler parameters has a physical inter- 
pretation. The Euler parameters q a correspond to an average orientation, in a geometric 
sense, of the two nodal cross-section orientations. This is demonstrated from the following 
example characterizing a state of constant curvature of a finite element shown in Figure 3. 
The orientation of the convected element frame is characterized by a rotation of an angle 
4> about an axis n a from the inertial reference frame, and the relative nodal cross-section 
orientations are characterized by a rotation from the convected frame of angles — r and r 
about axis nj, for nodes 1 and 2 respectively. The Euler parameters designating the total 
cross section orientation of the two nodes due to these combined effects can be expressed 

cos | cos ^ + n Q • n b sin | sin ^ 

— cos * sin ^ nj + cos j sin n a — sin f sin £ n a x n & 



9 2 


cos ^ sin j n j, 


cos 4 cos j — n a 


2 

+ 


cos j sin 4 


‘ (t> • T 

n& sin sin 
n a + sin ^ sin j n a 


x n& 


} 


(5.10) 


which is obtained by applying the quaternion product rule 44 to the Euler parameter defi- 
nitions 



of the relative nodal orientations and the convected orientation respectively. The average 
of the two nodal Euler parameters (5.10) is 



1 , 

2 ( 9i + ?2 ) 
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the norm of which is cos -j. When normalized, the above average is identical to the average 
orientation of the two nodes given by q a . It can be shown that for this example the dis- 
cretization (5.9) when substituted into (5.7) gives the finite element curvature computation 

4 . r 

* = J sin 2 "6 

which approximates the true curvature strain |rnft. The computation retains only the 
rotation parameters r originally defined relative to the rigid body orientation, and is thus 
invariant to the rigid body motions. For instances when the validity of the approximation 
is challenged, an incremental curvature computation can be made as discussed in the next 
section, from which the total curvature is obtained from an appropriate update procedure. 


5.2 Computation of the Strain Increments 


The strain increments axe defined from the virtual strains 
variational operator 8 with an incremental operator A as 


(3.12) by replacing the 


A7 


dAu 

~df 


+ 


0 

-A/?3 

A /?2 


1 


A/c 


d A/3 


such that A u and A/3 are finite analogs of the infinitesimal displacements and rotations 
8u and 8/3. For computation purposes, it becomes necessary to decompose the convected 
frame components of the virtual rotations of the of the cross-section 8/3 into a rotation due 
to rigid body motion 8q> and that due to deformations 8t as 


8/3 = 8r a + 8<f . 

This relation is derived by substituting the following definitions 

8/3 t = S T 8a T S , 8q t = SR R T 

8(p T = <5TT r , 8f T = 8SS t , 8tJ = S T 8f T S 
into the identity 


(5.11) 


R = S T 


6R = 8S T + S <!>T 


It is noted that 8a, 8ip, and St represent moving frame or spatial components referred to 
the defining reference frame, whereas 8/3 and 8r a represent material components referred 
back to the convected frame. From these definitions, the incremental strain A7 is given 
by 


A7 = T 


dAu 

~df 


+ 



+ 


0 

-At, 

At, 


<13 


) 


(5.12) 
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representing the membrane strain and transverse shear strain increments. Likewise, the 
incremental curvature representing the torsion and bending strains is given by 


A K 


d A r a 

dZ 


(5.13) 


as the incremental rotations A<p defined from the T matrix are constant over the element 
length. 


Essential for the use of these incremental strains is a proper definition and subsequent 
computation of the finite displacement and finite rotation increments. The incremental 
translations axe defined by 

A u = u n+1 — u" (5-14) 

as the displacements axe true vector quantities. The incremental rotations axe defined as 
follows. Rotations are updated by the product of orthogonal matrices via either 24 


R" +1 = R (/) R n = e 


R" 


8 t 

J 

n „© T 


= R n R( r ) = R n e 


(5.15) 


using the rotational vectors 9 or 0 based on the spatial or material reference frames 
respectively. It can be seen from the linearizations of the left and right rotational updates 24 

R n+1 ~ R n + $R - 

<5R = $ T R" = R" 0 r 


that the virtual rotations 


8<p T = 8T T t , 8tJ = S T 8S 

correspond to spatial and material rotation updates 

T n+1 = AT T n , S n+1 =. S"AS (5.16) 

respectively. Thus A (p and A r a axe defined as the rotational vectors parameterizing the 
matrices AT and AS respectively. Two different approximate methods which then extract 
this pseudovector from the given rotation matrix are used to obtain the incremental rota- 
tions. The particular approximation methods are chosen such that objective computations 
of the incremental strains (5.12) and (5.13) axe achieved. 

To this end, the first two terms of (5.12) 

-° 1 

-Av> 3 } (5.17) 
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must be computed such that the A <p rotation increment compensates for the rigid rotation 
contained in the displacement increment A u defined in (5.14). To accomplish this, A <p is 
computed by 

Aip T = AT" + 2 - AT n+ 2 T (5.18) 

where AT n+ 2 is defined from 

T n+ * = exp ( \a$ t ) T n = AT n+ 2 T n . (5.19) 

The computation was derived from the linear approximation 

T n+1 ~ ( I + A <p T ) T" 

rewritten as 

A£ r = ( T n+1 - T n ) T n+ * (5.20) 

and introducing (5.19) to achieve a skew symmetric matrix. In order to preserve rigid 
motions, the matrix T in the first term of (5.17) must be evaluated as T n+ 2. This is 
shown as follows from an example of the rigid rotation of an element in which A71 = 0. 
From (5.20), it is seen that the rotational term of (5.17) becomes 


( 0 1 

> = - T n+ 2 ( t” +1 - 



f Tn ) 


{ -A*>3 

- ) . 

u = 

{ Tl2 ( 

(5.21) 

l A<p 2 J 
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The finite element evaluation of the displacement term of (5.17) is given by 

<9Au 1 t x 

T-gr- = T J- ( Au 2 - Am ) (5.22) 

for the two-noded beam element of length £ e . For the rigid rotation of the second node 
about the first node, the incremental translational displacements are simply 

Auj = 0 , A u 2 = 4 ( ^ n+1 -tf) , ^ = ^ n+1 " 

as the direction cosines of the rotation are contained in the first row of the T matrix. 
Thus for (5.17) to be identically equal to zero, it is necessary to evaluate (5.22) using 
T n+ 2. To obtain the true stretch with respect to the neutral-axis reference frame at 
the current configuration, we simply rotate the mid-configuration computation up to the 
current configuration as 


AT n +2 


ipn+-^ 


dAu 

~w 
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/ ° 1 

< -Av?3 > 
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(5.23) 
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As in the preceding analysis, the incremental displacements for an arbitrary rotation and 
stretch axe given by 


Am = o , A u 2 = ( (£ e + d) t ” +1 - £ e t* ) 


where d represents a stretch relative to the original element length £ e . The rotational 
expression (5.21) remains valid, and the bracketed term in (5.23) becomes 


( - — ) T n+ J 

V L. + d ) 


i n+l 


£e + d 

Premultiplication of the above by AT n+ * results in the final computation 
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containing solely a measure of stretch regardless of the magnitude of the rigid rotation. 
The incremental rotations Ar a used to compute the remaining terms, 


A72 = 


0 

-Ar, 

At, 


a 3 


A/c = 


a 2 


d A r a 


(5.24) 


representing transverse shear and curvature strains respectively, are computed indepen- 
dently from Aip as follows. The rotation increments Ar a are obtained from the matrix AS 
defined in (5.16) denoting the relative orientation between the current deformation matrix 
S n+1 and the past deformation matrix S n rigidly rotated to the current convected frame. 
Another method to extract a rotation pseudovector from a given orthogonal rotation ma- 
trix given by 43 

X~T = z ( AS, - AS l ) 
a •' 1 + tr AS, 


1 = 1,2 


(5.25) 


is used to define Ar a at each element node. The above method yields a simpler and 
more accurate computation of a rotation vector than (5. IS). Whereas (5. IS) was necessary 
to compute A<p such that the rigid rotations within (5.17) are preserved, (5.25) is used 
within (5.24) as this computation is made from matrices which by construction contain 
information solely due to deformation. Given the nodal rotation increments, the locking- 
free form of the elemental shear strain is obtained from the nodal average as 
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and the elemental curvature is computed from the finite element approximation 

Ak = j (A r a2 - Ar 0l ) 

This completes the computational procedures for the incremental strains. The detailed 
strain computations of (5.12) and (5.13) are used in (5.4) to determine the stress incre- 
ments, from which the current stress state is obtained from the update procedure (5.3). 


6. Numerical Examples 

The computational techniques, namely the staggered multibody dynamics solution 
procedure combining the generalized coordinate integrator and the constraint force solver 
discussed in Section 4 and the finite element computations of the beam internal force 
discussed in Section 5, have been implemented into a Fortran 77 software package. The 
result is a robust method which solves the present formulation of the equations of motion 
of an arbitrary assemblage of flexible beams and rigid bodies. In order to demonstrate the 
current software capabilities, the following examples highlighting the flexible motion of the 
beam component axe presented. 

The first example is included to verify the geometric stiffening phenomena exhibited 
by a rotating beam 6,18,21,28 . The beam is pinned at the left end; the other end remains 
free. The following material and geometric properties were used: 

EA = 2.8 x 10 7 lb, GA = 1.0 x 10 7 lb, El = 1.4 x 10 4 76 in 2 

pA = 1.2 Ibm/in , pi = 6.0 x 10 -4 Ibm in, l = 10 in. 

A prescribed angular rotation about the e 3 axis of 

e(t ) = I + (cos W " 1)1 md 0 £ ‘ - 15 

\ ( 6t — 45) rad t > 15 sec 

is applied at the pinned end. The time history of the tip deflection relative to a refer- 
ence frame coinciding with the prescribed angular position and the time history of several 
configurations of the beam are given in Figure 4. As alluded to in the introduction, an 
overall steady rotation of the beam gives rise to a centrifugal force which is responsible 
for a change in the bending stiffness that cannot be predicted using linear deformation 
theories. After initial increasing tip deflections, the beam begins to stiffen as the angular 
velocity increases due to the centrifugal inertia force. As the angular velocity reaches a 
constant state, the beam then reaches a steady state phase of small vibrations. This ex- 
ample shows the capability of the nonlinear strain formulation to automatically account 
for the geometric stiffening effect. The results are comparable to those presented by Simo 
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and Vu-Quoc 28 . To reproduce these results with alternative methods as the substructur- 
ing technique 21 , a convergence analysis based on the selection of mode shapes must be 
performed. 

The next examples exhibit the combined large deformation and large rotation capa- 
bilities of the present formulation. In the first instance, the beam is pinned as above and 
is subjected to given initial velocity impulses exciting various deformation mode shapes 
under planar motion. The following material and geometric properties were used in order 
to witness finite deformations: 


EA = 4.0 x 10 7 lb, GA = 2.0 x 10 7 lb, El = 1.3 x 10 7 lb in 2 

pA = .98 Ibm/in, pi = 3.3 x 10 -2 Ibm in, l = 200 in. 

The initial velocity profiles with the resulting time histories of several deformed config- 
urations are given in Figures 5, 6, and 7. It is noted the versatility of the formulation 
in its ability to capture the response to a variety of situations in which different funda- 
mental modes of the beam are excited. The approach avoids the difficulty of tailoring 
the selection of modes shapes of the flexible components to the given problem at hand. 
The repeatability of the deformation shapes through time is due to the invariance of the 
internal force computations to the overall rigid motion. This property of computational 
objectivity is further illustrated in Figure 8 which shows the time history of the strain, 
kinetic, and total energy over four revolutions for the first bending mode example. The 
nature of the time integration and internal force algorithms are such that the conservation 
of energy is retained computationally, as seen by the fact that the total energy remains 
constant over all the revolutions. Similar results, not presented within, are obtained for 
the other deformation examples. 


To present the applicability of the flexible beam component within the multibody 
dynamics framework, the final example of a spatial double pendulum is given. The double 
pendulum is modeled with two beams; a spherical joint connects the last node of the first 
beam to the first node of the second beam and also pins the first node of the first beam. 
It is noted that the joint connection can easily be accounted for from a finite element 
assemblage which leaves the rotational degrees of freedom free at the hinge location. The 
method was used to verify the results obtained using the Lagrange multiplier solver on the 
augmented equations described in Section 3.5. In the first case, the pendulum is subjected 
to a gravity field in the vertical z-direction and an initial velocity impulse in the horizontal 
x-y plane such that soley rigid motion is excited. The problem is run for four cases of 
increasing beam flexibility as follows: 


1. EA = 1.0 x 10 4 lb 

2. EA = 1.0 x 10 3 lb 

3. EA = 2.0 x 10 2 lb 

4. EA = 1.0 x 10 2 lb 

with the remaining parameters 


GA = 0.5 x 10 4 lb 
GA = 0.5 x 10 3 lb 
GA = 1.0 x 10 2 lb 
GA = 0.5 x 10 2 lb 
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pA = 1 Ibm/in pi = .833 x 10 -3 Ibm in 1 = 1 in 
held constant. The initial velocity impulse, and the spatial trajectories of the mass 
center of the second beam as projected on the x-y and x-z planes is given in Figure 9. The 
trajectory of the first case coincides exactly with a rigid body solution to the problem, 
and the slight deviation of the trajectories due to the increasing flexibility can be seen. 
The energy time histories for the problem, given in Figure 10, verify the computational 
objectivity of the algorithm as again energy is identically conserved. Again, the invariance 
of the internal force calculations in the three dimensional environment is witnessed by the 
negligible strain energy contribution for all of the flexible cases. The time integration of 
the spatial kinematics preserves the balance between the kinetic and potential energies of 
the problem. Next, the flexible double pendulum is given an initial velocity impulse to 
excite deformation motion as well as the rigid motion. For this case the parameters used 
were 

EA = 1.8 x 10 7 8 lb, GA = 0.9 x 10 8 lb, El = 1.4 x 10 8 lb in 2 

pA = .98lbm/in, pi = 0.67 Ibm in, l = 120 in. 

The initial velocity profile, the resulting time histories of several deformed configurations 
and energy time history are given in Figure 11, exhibiting the large spatial rotation and 
deformation capabilities of the formulation. The energy conservation is retained for the 
computations of spatial deformations. 

Further examples of large scale multibody systems are in process, and these results 
are to be presented in the near future. 


7. Concluding Remarks 

A flexible beam finite element that is readily incorporated into multibody dynamics 
applications has been presented. The beam formulation is based on fully nonlinear strain 
measures which remain invariant to rigid body motions. The model retains a Cauchy 
stress and physical strain description, and as such it can be easily interfaced with real- 
time slewing control applications as the measured strains can directly be used as a feedback 
signal without requiring sophisticated transformations. In addition, the formulation uses 
an inertial reference for the beam dynamics such that the degrees of freedom of the flexible 
component are defined in the same sense as the rigid components by including without dis- 
tinction both the rigid and flexible deformation motions. The consequence is adaptability 
into multibody dynamics methodologies as numerical solution procedures for the integra- 
tion of spatial kinematic systems can directly be applied to the generalized coordinates of 
both the rigid and flexible components. The success of the approach relies on an accurate, 
computation of the nonlinear internal force term. For this reason, the calculation of finite 
strain increments has been presented which are invariant to arbitrary rigid motions of the 
beam. The proposed methodology is suitable to treat the dynamics of flexible beams which 
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undergo a variety of structural deformations in addition to the large overall motions. The 
same approach can be used in formulating other types of structural components. 
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Figure 4. Geometric Stiffening ( 

(a) Tip Deflection vs. Time 

(b) Displacement History 
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Figure 5. First Bending Mode ( 8 Elements ): 

(a) Initial Beam Position vs. Initial Velocity Profile 

(b) Displacement History 
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Figure 6. Second Bending Mode ( 12 Elements ): 

(a) Initial Beam Position vs. Initial Velocity Profile 

(b) Displacement History 









Figure 9. Spatial Double Pendulum ( 16 Elements ): 
(a) Second Beam Trajectory: X-Y Plane 

. (b) Second Beam Trajectory: X-Z Plane 
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